******************
*** 1_PART ONE ***
*** Simulation ***
******************
cd "~"
frame create simu
frame change simu

clear all
set more off
version 19.5
set seed 13579

* Input predicted constraint levels (see file "2a_Part I")
input str6 system str12 period byte K double rho_start double rho_end double se_rho_start double se_rho_end
"BSA"   "22 vs. 02" 4   0.0875048   0.1412387   0.0038848   0.0037695
"BSA"   "12 vs. 02" 4   0.0977978   0.1037745   0.0040746   0.0037507
"BSA"   "22 vs. 12" 4   0.0993402   0.1536674   0.0045535   0.0045375
"BSA"   "22 vs. 16" 4   0.1237224   0.151166    0.0066061   0.0057223

"BSec"  "22 vs. 02" 2   0.1365416   0.1598804   0.0076763   0.0074506
"BSec"  "12 vs. 02" 2   0.1422674   0.1383483   0.0079519   0.0073246
"BSec"  "22 vs. 12" 2   0.1322014   0.1723463   0.0094728   0.0094359
"BSec"  "22 vs. 16" 2   0.1525602   0.1718867   0.0130016   0.0111012

"BSB"   "22 vs. 16" 5   0.1349447   0.1853219   0.0081527   0.0073932
end

* Settings
local N      = 20000
local anchor = 10

* Prep result columns in the MAIN frame
gen double p_start = .
gen double p_end   = .
gen double ratio   = .
gen double est_end = .

quietly count
local R = r(N)

* Create a SEPARATE frame for simulations
capture frame drop __sim__
frame create __sim__

forvalues i = 1/`R' {

    local Ki   = K[`i']
    local rhos = rho_start[`i']
    local rhoe = rho_end[`i']

    * ------- START scenario in sim frame -------
    frame change __sim__
    clear
    matrix Rm = J(`Ki',`Ki',`rhos')
    forvalues j = 1/`Ki' {
        matrix Rm[`j',`j'] = 1
    }
    * explicit varlist z1 z2 ... zK
    local varlist
    forvalues j = 1/`Ki' {
        local varlist `varlist' z`j'
    }
    quietly corr2data `varlist', n(`N') corr(Rm)
    gen byte all_pos=1
    gen byte all_nonpos=1
    forvalues j = 1/`Ki' {
        replace all_pos    = all_pos    & (z`j' >  0)
        replace all_nonpos = all_nonpos & (z`j' <= 0)
    }
    gen byte fully = (all_pos | all_nonpos)
    quietly summarize fully
    local ps = r(mean)

    * ------- END scenario in sim frame -------
    clear
    matrix Rm = J(`Ki',`Ki',`rhoe')
    forvalues j = 1/`Ki' {
        matrix Rm[`j',`j'] = 1
    }
    * explicit varlist again
    local varlist
    forvalues j = 1/`Ki' {
        local varlist `varlist' z`j'
    }
    quietly corr2data `varlist', n(`N') corr(Rm)
    gen byte all_pos=1
    gen byte all_nonpos=1
    forvalues j = 1/`Ki' {
        replace all_pos    = all_pos    & (z`j' >  0)
        replace all_nonpos = all_nonpos & (z`j' <= 0)
    }
    gen byte fully = (all_pos | all_nonpos)
    quietly summarize fully
    local pe = r(mean)

    * ------- Write results back in MAIN frame -------
    frame change default
    replace p_start = `ps' in `i'
    replace p_end   = `pe' in `i'
    replace ratio   = p_end / p_start in `i'
    replace est_end = 10 * ratio in `i'
}

* Done with sim frame
frame drop __sim__

*------------------------------------------------------
* Contrast-specific baselines (baseline = row's start)
*------------------------------------------------------
drop ratio est_end
gen double p_base  = p_start
gen double ratio   = p_end / p_base
gen double est_end = 10 * ratio
gen double diff    = est_end - 10

label var p_base  "Contrast baseline p (row-specific)"
label var ratio   "End/Baseline ratio (row-specific)"
label var est_end "Anchored end per 100 (baseline=10)"
label var diff    "Increase vs baseline (per 100)"
format est_end diff %6.2f

*------------------------------------------------------------------
* Delta-method style 95% CIs for est_end and diff
*------------------------------------------------------------------

* Row-specific baseline rho and SE (just take rho_start and se_rho_start)
gen double rho_base    = rho_start
gen double se_rho_base = se_rho_start

* storage
gen double se_est_end = .
gen double est_end_lo = .
gen double est_end_hi = .
gen double diff_lo    = .
gen double diff_hi    = .

* small step for numerical derivative dp/drho
local h = 0.001

* ---- helper program (define ONCE) ----
capture program drop __getp
program define __getp, rclass
    version 19.5
    // returns r(p) = Pr(all-consistent) for K items under equicorr(rho)
    args K rho N
    tempname Rm
    matrix `Rm' = J(`K',`K',`rho')
    forvalues j = 1/`K' {
        matrix `Rm'[`j',`j'] = 1
    }
    // explicit varlist: z1 z2 ... zK
    local vlist
    forvalues j = 1/`K' {
        local vlist `vlist' z`j'
    }
    quietly corr2data `vlist', n(`N') corr(`Rm')
    gen byte all_pos=1
    gen byte all_nonpos=1
    forvalues j = 1/`K' {
        replace all_pos    = all_pos    & (z`j'>0)
        replace all_nonpos = all_nonpos & (z`j'<=0)
    }
    gen byte fully = all_pos | all_nonpos
    quietly summarize fully
    return scalar p = r(mean)
    drop _all
end

* temp frame for quick evaluations
capture frame drop __deriv__
frame create __deriv__

quietly count
local R = r(N)

forvalues i = 1/`R' {
    * pull row values (as locals)
    quietly su K in `i', meanonly
    local Kvar = r(mean)

    quietly su rho_end in `i', meanonly
    local re_hat = r(mean)
    quietly su se_rho_end in `i', meanonly
    local se_re = r(mean)

    quietly su rho_base in `i', meanonly
    local rb_hat = r(mean)
    quietly su se_rho_base in `i', meanonly
    local se_rb = r(mean)

    * clamp rhos to [0,.95] and build +h steps
    if (`re_hat' < 0)   local re_hat = 0
    if (`re_hat' > .95) local re_hat = .95
    if (`rb_hat' < 0)   local rb_hat = 0
    if (`rb_hat' > .95) local rb_hat = .95
    local reh = `re_hat' + `h'
    if (`reh' > .95)    local reh = .95
    local rbh = `rb_hat' + `h'
    if (`rbh' > .95)    local rbh = .95

    * evaluate p() at end and end+h
    frame change __deriv__
    clear
    __getp `Kvar' `re_hat' 20000
    local p_end = r(p)
    clear
    __getp `Kvar' `reh' 20000
    local p_end_h = r(p)

    * evaluate p() at base and base+h
    clear
    __getp `Kvar' `rb_hat' 20000
    local p_base = r(p)
    clear
    __getp `Kvar' `rbh' 20000
    local p_base_h = r(p)

    frame change default

    * numerical derivatives dp/drho
    local dpdr_end  = (`p_end_h'  - `p_end')  / `h'
    local dpdr_base = (`p_base_h' - `p_base') / `h'

    * delta-method SEs for p_end and p_base
    local se_p_end  = abs(`dpdr_end')  * `se_re'
    local se_p_base = abs(`dpdr_base') * `se_rb'

    * ratio r = p_end / p_base ; var(r) ≈ (se_p_end/p_base)^2 + (p_end*se_p_base/p_base^2)^2
    local r_hat = `p_end' / `p_base'
    local se_r  = sqrt( ( `se_p_end' / `p_base' )^2 + ( (`p_end' * `se_p_base') / (`p_base'^2) )^2 )

    * est_end and CI
    local se_est = 10 * `se_r'
    replace se_est_end = `se_est' in `i'
    replace est_end_lo = est_end[`i'] - 1.96*`se_est' in `i'
    replace est_end_hi = est_end[`i'] + 1.96*`se_est' in `i'

    * diff = est_end - 10 ; CI uses same SE
    replace diff_lo = diff[`i'] - 1.96*`se_est' in `i'
    replace diff_hi = diff[`i'] + 1.96*`se_est' in `i'
}

frame drop __deriv__

label var est_end_lo "est_end 95% CI lo"
label var est_end_hi "est_end 95% CI hi"
label var diff_lo    "diff 95% CI lo"
label var diff_hi    "diff 95% CI hi"

encode system, gen(sys1)
recode sys1 (1=1 "Belief System {it:A}") (2=3 "Belief System {it:B}") (3=2 "Economic Belief System"), gen(bsys)

encode period, gen(per1)
recode per1 (1=1 "{it:Ref.}: 2002") (2=.) (3=2 "{it:Ref.}: 2012") (4=3 "{it:Ref.}: 2016"), gen(per)

** Figure 3b
twoway (bar diff bsys, by(per, rows(1) note(" ")) bc(edkblue) fcolor(edkblue%0)) (rcap diff_lo diff_hi bsys, lcolor(edkblue*.5) lwidth(thin)), ///
	subtitle(,bcolor(white)) ///
	xlab(.5 " " 1 "BS {it:A}" 2 "Ec. BS" 3 "BS {it:B}" 3.5 " ", nogrid) xtitle(" ") ///
	ylab(-.5 " " 0 "0" .5 " " 1 "1" 1.5 " " 2 "2" 2.5 " " 3 "3", angle(0)) ///
	yline(0, lc(cranberry*.5) lwidth(vthin) lp(shortdash)) ///
	legend(rows(1) order(1 "Increase vs baseline" 2 "95% CI") pos(6) region(lstyle(none))) scheme(s1mono)
** remove ticks in gr editor